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ABSTRACT 


A  proposal  has  been  made  to  compute  a  2-D  numerical  solution  to 
the  stream  function  -  vorticity  equation  to  analyze  the  flow  field,  both 
of  an  airfoil  in  ground  effect  and  of  an  air  cushion  vehicle  of  arbitrary 
geometry  supported  by  peripheral  jets.  This  report  is  an  investigation  of 
both  problems.  Finite  difference  solutions  of  the  stream  function-vorti- 
city  equation  are  computed  on  a  body  fitted  coordinate  system.  The 
Laplace  equation  for  streamlines  around  an  airfoil  in  ground  effect  is 
solved,  and  results  are  compared  to  an  analytic  finite  singularity  solu¬ 
tion.  Comparison  of  C^  and  C^  at  varying  angles  of  attack  and  ground 
distance  indicate  good  agreement.  The  finite  difference  method  is  found 
to  be  sensitive  to  grid  point  distribution  about  the  leading  edge,  for 
which  a  criterion  is  established.  Since  an  air  cushion  vehicle  acts  as  a 
flow  source  within  the  field,  the  finite  difference  expression  is  modified 
to  accommodate  a  branch  cut.  The  presence  of  vorticity  is  essential  to 
definition  of  the  peripheral  jets.  The  vorticity  distribution,  derived 
from  arbitrary  velocity  distribution  at  jet  outlets  and  assumed  to  be 
constant  along  streamlines,  leads  to  a  stable  solution  of  the  stream 
function-vorticity  equation. 
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2-D  FLOW  NUMERICAL  SOLUTION  FOR 


AIRFOIL  AND  HOVERCRAFT  IN  GROUND  EFFECT 

I.  INTRODUCTION 


Background 

The  Air  Force  Flight  Dynamics  Laboratory  is  currently  investigating 
the  feasibility  of  designing  a  high  performance  aircraft  with  air  cush¬ 
ioned  wings  instead  of  standard  landing  gear  to  achieve  a  take-off  and 
landing  capability  on  rough  and  oomb-damaged  runways.  The  wing  air  cush¬ 
ion  is  to  be  generated  by  a  peripheral  jet.  Investigating  the  behavior  of 
such  a  design  is  a  complicated  fluid  mechanics  problem  that  combines  air 
cushion  dynamics  with  a  wing  in  ground  effect.  Current  methods  for  air¬ 
foil  flow  field  solution  range  from  simple  potential  flow  solutions  to 
viscous,  compressible  Navier-Stokes  solutions.  Air  cushion  flow  fields 
are  determined  through  experimental  techniques.  Current  air  cushion  anal¬ 
ysis  (Ref  3,  7)  is  based  upon  assuming  a  simple  jet  geometry  to  enable 
integration  of  pressure  gradients.  In  the  commonly  used  exponential  theory 
jet  streamlines  with  a  constant  radius  of  curvature  are  assumed,  and  pres¬ 
sure  gradients  across  streamlines  supply  the  centripital  acceleration. 

These  methods  provide  reasonable  results  for  symmetric  cases  only. 

Improved  flow  analysis  of  air  cushion  dynamics  would  diminish  the  need  to 
evaluate  design  parameters  exclusively  through  expensive  experiments.  A 
proposed  approach  is  to  determine  the  flow  field  for  an  air  cushioned  air¬ 
foil  by  using  2-D  finite  difference  solution  to  the  stream  function- 
vorticity  equations. 
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Purpose  of  This  Study 


The  general  purpose  is  to  obtain  a  generally  applicable  method  for 
solution  of  air  cushioned  wing  flow  field.  In  particular,  a  2-D  flow 
field  solution  by  using  the  stream  function-vorticity  equations  is  pro¬ 
posed.  This  study  is  a  small  scaled  feasibility  study  to  the  proposed 
solution  and  includes  two  separate  aspects  of  the  problem.  The  first  as¬ 
pect  is  investigating  the  behavior  of  a  finite  difference  numerical  solu¬ 
tion  for  an  airfoil  in  ground  effect  and  comparing  the  results  to  those 
obtained  by  an  independent  method.  The  second  aspect  is  to  solve  a  sym¬ 
metric  case  of  air  cushion  vehicle  with  peripheral  jet  and  obtain  a 
physically  reasonable  flow  field. 

Procedure 

As  mentioned  the  flow  fields  in  this  study  are  solved  by  using  the 
stream  function-vorticity  formulation.  Use  of  the  stream  function  guar¬ 
antees  conservation  of  mass.  The  general  incompressible  vorticity  equation 
that  is  derived  from  momentum  equations  has  the  form 

=  f^  +  ^  *  7(0  =  v  •  (1) 

Equation  (1)  indicates  that  in  the  absence  of  diffusion  the  vorticity 
must  be  constant  along  streamlines,  in  which  case  equation  (1)  can  be  re¬ 
written  in  stream  function  formulation  for  planar  flow  as 

v2t|>  =  -co(i p)  (2) 

A  tacit  assumption  throughout  the  study  is  that  neither  the  initial  dis¬ 
tribution  of  vorticity  in  the  peripheral  jet  nor  the  vorticity  transport 
has  a  critical  impact  on  pressure  distribution  under  the  vehicle  (both 
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points  must  be  checked  by  an  additional  feasibility  study).  In  summary, 
the  vorticity  is  determined  from  an  arbitrary  assumed  velocity  distri¬ 
bution  across  the  peripheral  jet  outlets  in  the  air  cushion  case.  For 
the  case  of  the  airfoil  in  ground  effect  the  vorticity  is  set  to  zero. 

The  numerical  formulation  consists  of  generation  of  body  fitted 
coordinate  system  by  solving  Poisson  equations  for  transformed  coordi¬ 
nates  (C,n)  as  a  function  of  physical  coordinates  (x,y).  Free  parameters 
in  grid  generation  equations  as  well  as  freedom  in  the  selection  of  boun¬ 
dary  points  location  allow  the  manipulation  of  the  physical  grid  to  a 
configuration  which  provides  enhanced  resolution  where  it  is  needed.  The 
stream  function-vorticity  equations  and  the  body  fitted  coordinate  system 
equations  are  both  transformed  to  the  U»n)  plane  and  solved  by  SOR  finite 
difference  techniques. 

Specific  Problems  Addressed 

1.  Cylinder  in  uniform  potential  flow 

This  case  is  solved  to  verify  both  the  body  fitted  coordinate 
system  generation  and  the  potential  flow  stream  function  solu¬ 
tion  by  comparison  with  well  known  analytical  results. 

2.  Airfoil  in  free  stream  and  ground  effect 

The  2-D  flow  field  numerical  solution  is  used  even  though  it  is 
realized  that  some  important  aspects  of  wing  in  ground  effect 
are  related  to  3-D  analysis,  such  as  the  dramatic  decrease  in 
induced  drag  and  large  differences  in  handling  qualities.  The 
airfoil  flow  field  is  analyzed  with  both  a  finite  difference 
method  and  a  finite  singularity  method.  Lift  and  quarter  chord 
moment  coefficients  variation  with  angle  of  attack  and  ground 
distance  are  derived  by  both  methods  and  compared. 
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3.  Air  cushion  device  with  peripheral  jets 

Introduction  of  branch  cut  to  relieve  multi-valued  stream 
function  from  fluid  source  within  the  vehicle.  To  obtain  a 
well  defined  jet,  vorticity  is  introduced  into  shear  layers 
at  the  edges  of  the  peripheral  jets  and  assumed  to  remain 
constant  along  each  streamline. 
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II.  FINITE  SINGULARITY  METHOD 


In  this  study,  the  airfoil  in  ground  effect  is  solved  for  an  in- 
vicid  incompressible  flow  using  a  finite  difference  method.  This  chapter 
describes  an  airfoil  in  ground  effect  solution  using  a  finite  singularity 
method  which  establishes  a  comparison  standard  for  the  finite  difference 
solution. 

Thin  Airfoil  Theory 

The  problem  of  calculating  the  flow  field  and  the  aerodynamic  prop¬ 
erties  of  an  arbitrary  airfoil  with  no  restriction  as  to  its  thickness, 
camber,  or  angle  of  attack,  is  complex  in  practice.  Thin  airfoil  theory 
provides  simplification  of  the  mathematical  problem  with  a  rather  wide 
applicability.  In  thin  airfoil  theory  the  flow  around  the  airfoil  is 
represented  by  a  sheet  of  vortices  lying  along  the  camber  line.  The  per¬ 
turbed  flow  field  is  determined  by  the  following  (Ref  1): 

1.  Differential  equation 


2.  Boundary  conditions 

(U  +  Vip)  ♦  vF  =  0  on  airfoil  (4) 

v<j>  «♦.  o  at  infinity  (5) 

3.  Kutta  condition 

The  velocity  is  finite  and  continuous  at  the  trailing  edge. 

Once  the  airfoil  shape  is  replaced  by  a  sheet  of  vorticity  of 
strength  y(x)  at  the  camber  line  and  thin  airfoil  assumptions  are  applied, 
the  mathematical  problem  reduces  to  one  singular  integral  equation. 


5 


(6) 


dYa  _  1 

dx  2n 


Cl  2 

-l 


llii 

x-C 


d£ 


Solution  of  the  integral  equation,  to  determine  y(x),  can  be  done 
by  the  following  techniques: 

1.  Sohngen  integral  inversion 

2.  Prandtl  Glauert  transformation 

3.  Finite  singularity  method 

Given  the  vorticity  distribution  y  along  the  camber  line  (or  along 
the  x-axis),  the  lift  and  moment  coefficients  are  easily  found  by  inte¬ 
gration. 


C 


i 


C/2 

/  y(x)dx 
-c/2 


2  w  / 

w-rz-  /  (x.-x)  '  Y(x)dx 

Voc  -C/2  0 


(7) 

(8) 


The  rest  of  this  chapter  deals  with  finite  singularity  solution  for 
a  flat  plate  in  ground  effect.  The  solution  provides  y(x),  C£,  Cm  as 
function  of  ground  distance  and  angle  of  attack. 

Method  Description 

The  finite  singularity  method  is  an  approximate  numerical  method 
to  solve  the  integral  equation  (6).  A  finite  number  of  point  vortices 
is  obtained  by  solving  a  set  of  algebraic  equations,  which  is  a  very 
convenient  computational  formulation  (Ref  2). 

The  airfoil  is  divided  into  N  elements  (Fig  1).  Increasing  N  im¬ 
proves  the  accuracy,  but  usually  ten  elements  are  enough.  A  vortex  of 
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finite  strength  is  placed  at  the  25%  of  each  element,  points  (X^ , ) . 
The  tangential  flow  boundary  condition  is  satisfied  at  75%  of  each  ele¬ 
ment,  points  (X - , Y . ) .  The  set  of  algebraic  equations  to  be  solved  are: 

J  J 


u  fa 

03X 


Xj’yj 


•1  £  WYiXXj-V 

=  W(Xj,yj)  -  -  2^  i=l  rXj-X1^+(Yj-Y1):f  0  «  1,2,. 


Flat  Plate  in  Ground  Effect 

As  mentioned  above,  this  method  is  to  be  applied  to  a  thin  symmetric 
airfoil  in  ground  effect  for  the  purpose  of  establishing  a  point  of  com¬ 
parison  for  verification  of  the  finite  difference  code.  In  the  context 
of  thin  airfoil  theory  a  1%  thickness  symmetric  airfoil  appears  as  a  flat 
plate.  Ground  effect  in  the  context  of  the  finite  singularity  method  is 
established  by  using  an  image  flat  plate  to  satisfy  tangential  flow  on  the 
ground  (Fig  2).  The  tangential  flow  condition  at  point  (x.,y.)  in  the 

J  J 

plate  is  determined  by  vortices  on  both  the  plate  and  the  image,  which 
leads  to  the  following  equation. 


u|ya 

03X 


i  JL  ri,(xj-xi) 

2,  £  T 


yryi 


1  jt  ri'(xj'xi) 

!+  57  Si  <xrv"+lV,i 


Substituting  the  following  relations  into  equation  (10) 


X  -i  *cosa  <  =  i  or  j 

K  1C 

Y 

K 


(ID 

(12) 


a-i  COSa 

K 


tc  =  i  or  j 


(13) 


gives: 


-2nU0a 


r^'COSa  £  ■  I\(;Aj.-A..)c0Sa 

.  ~  ”  i  =1  a7-a7  +  1=1  (a  --a.  )^cos2a+[2a-(A  .+jt.)sina;F 

J  1  J  1  J  • 


04) 


Where  points  a., a.  are  calculated  for  N  equal  elements  as: 

■  J 

Ai  =  (-0.5  +  i^I+  ^-)*  C  i  =  1 . N 

l.  =  (-0.5  +  2^).  C  j  =  1 . N 

Equation  (14)  is  written  in  a  matrix  form 


05) 

06) 


-2ttUoci  =  COSa  •  [A]{r.j} 


(17) 


The  vortex  strengths  are  obtained  by  inverting  matrix  A. 
fr  \  =  "foPoPf.  fA”1! 

Ur  COSa  LM  J 


(18) 


Integrating  equations  (7)  and  (8)  numerically  provides  the  lift  and  mom¬ 
ent  coefficients 

N 


2*COSa  *- 

Vc  fl  1 


(19) 


r  _  2c0Sa  /C  „  ^  „ 

cm  [j  T7  (t  -*■,)•  r.- 

c/4  uou  i=l  *  1  1 


(20) 


Equations  (18),  (19),  and  (20)  are  solved  numerically  for  20  equal  length 
elements. 

Results 

Tabulated  and  graphic  representation  (Tables  1  and  2,  Figs  3r6)  of 
lift  and  quarter  chord  moment  coefficients  as  function  of  ground  distance 
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and  angle  of  attack,  are  shown.  The  results  indicate  an  increase  in 
lift  coefficient  and  pitch  down  quarter  chord  moment  coefficient,  for  a 
given  angle  of  attack  as  the  airfoil  approaches  the  ground.  The  increase 
in  lift  and  pitch  down  quarter  chord  moment  coefficients  start  to  be 
significant  when  the  ground  distance  equals  the  chord  length  and  less. 

The  obtained  results  will  be  compared  against  those  obtained  by  the 
finite  difference  method  in  Chapter  V. 


GROUND-O I  ST flNCE/CHORO 


FIGURE  4.  Lift  Coefficient  vs.  Ground  Distance 
Angle  of  Attack  as  a  Parameter. 


A  -  10°,  £  -  4° 


ANGLE  OF  ATTACK 


FIGURE  5.  Quarter-Chord  Moment  Coefficient  vs.  Angle  of  Attack 
Ground  Distance  as  a  parameter. 

0  -  0.25c  A  -  0.50c 


+  -  1 .  OOr. 


X  -  Free  stream 


TABLE  1.  Lift  Coefficient  as  Function  of  Angle  of  Attack 
and  Ground  Distance  -  Finite  Singularity  Method 


a 

GROUND  DISTANCE 

INFINITY 

1.00 

0.50 

0.25 

0.0 

0.00 

0.00 

0.00 

0.00 

1.0 

.11 

.12 

.12 

.17 

2.0 

.22 

.23 

.25 

.33 

3.0 

.33 

.35 

.39 

.50 

4.0 

.44 

.46 

.52 

.67 

5.0 

.55 

.58 

.65 

.84 

6.0 

.66 

.70 

.78 

1.01 

7.0 

.77 

.81 

.92 

1.18 

8.0 

.88 

.93 

1.05 

1.35 

9.0 

.99 

1.04 

1.18 

1.52 

10.0 

1.10 

1.16 

1.31 

1.70 

TABLE  2.  Quarter  Chord  Moment  Coefficient  as  Function 
of  Angle  of  Attack  and  Ground  Distance  - 
Finite  Singularity  Method 


GROUND  DISTANCE 

a 

INFINITY 

1.00 

0.50 

0.25 

0.0 

0.00000 

0.00000 

0.00000 

0.00000 

1.0 

-.00000 

-.00076 

-.00237 

-.00578 

2.0 

-.00000 

-.00151 

-.00470 

-.01129 

3.0 

-.00000 

-.00226 

-.00698 

-.01654 

4.0 

-.00000 

-.00300 

-.00922 

-.02154 

5.0 

-.00000 

-.00374 

-.00142 

-.02629 

6.0 

-.00000 

-.00446 

-.01357 

-.03079 

7.0 

-.00000 

-.00519 

-.01 563 

-.03504 

8.0 

-.00000 

-.00590 

-.01777 

-.03905 

9.0 

-.00000 

-.00662 

-.01982 

-.04279 

10.0 

-.00000 

-.00732 

-.02183 

-.04627 
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III.  COORDINATE  TRANSFORMATION 


Coordinate  transformation  (or  generation  of  boundary  fitted  coordi¬ 
nate  system)  is  a  general  method  that  can  be  used  for  multiple  bodies  with 
arbitrary  boundaries.  The  method  can  be  used  for  solution  of  any  type  of 
partial  differential  equations.  The  mathematical  background  and  develop¬ 
ment  are  described  in  detail  in  Ref  4.  In  this  work  the  method  is  used  to 
solve  the  stream  function  equation  for  various  boundaries  and  boundary  con¬ 
ditions.  In  this  chapter  a  very  brief  discussion  of  the  method  is  introduced. 
Mathematical  Development 

For  purposes  of  discussion  in  this  work  the  region  of  interest  is 
two  dimensional  and  doubly  connected.  The  solution  of  all  the  partial  dif¬ 
ferential  equations  is  done  in  the  transformed  plane  which  is  rectangular 
in  shape  and  has  a  square  mesh.  Figure  7  illustrates  an  arbitrary  two  di¬ 
mensional  doubly  connected  region  with  arbitrary  inner  and  outer  boundaries 
that  is  transferred  into  a  rectangular  region.  The  inner  contour  maps  into 
the  ni  line  and  the  outer  contour  maps  into  the  n2  line.  The  contours 
connecting  the  inner  and  outer  boundaries  map  into  £  .  and  £  .  Since 

they  coincide  in  the  physical  plane,  they  constitute  re-entrant  boundaries 
in  the  transformed  plane.  This  type  of  mapping  is  used  for  all  transfor¬ 
mations  in  this  work.  The  coordinates  of  each  grid  point  in  the  physical 
plane  are  generated  as  solutions  of  Poisson  equations  in  the  physical  plane, 
with  Dirichlet  boundary  conditions  on  all  boundaries.  Control  functions  P 
and  Q  appear  as  inhomogeneous  terms  in  the  Poisson  equations,  that  are: 


Physical  Plane 
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Transformed  Plane 


Figure  7.  Field  Transformation  -  Single  Body 
(from  Ref  4,  p.  170) 
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«xx*eyy’  p(5-">  (21> 

Qxx  +  5yy  =  (22) 

This  system  in  the  transformed  form  becomes: 

c»Xrr  -  2gXr  +  YX  +  J2  (PXr  +  QX  )  =  0  (23) 

44  4n  nn  4  n 

aY  -  2gY  +  yY  +  J2  (PYc  +  QY  )  =  0  (24) 

44  4n  nn  4  n 


Where  a,  g,  y,  J  are  functions  of  first  order  partial  derivatives 
of  the  physical  coordinates  (X,Y)  with  respect  to  transformed  coordinates 
(5,n).  This  type  of  mapping  physical  boundaries  into  the  transformed  plane 
is  arbitrary  and  not  limited  to  the  mapping  described  above. 

Coo r d i na t e  System  Con troJ_ 

Spacing  of  coordinate  lines  on  the  boundaries  is  accomplished  easily 
as  boundaries  points  are  given  as  input  data.  Also,  the  coordinate  lines 
that  map  into  4max  and  cmin  in  the  transformed  plane  are  easily  controlled 
by  input  data.  However,  a  spacing  of  coordinate  lines  within  the  field  is 
controlled  by  varying  the  elliptic  generating  system  (changing  P  and  Q). 

In  the  control  system  used,  P  is  zero  and  Q  is  a  sum  of  exponential 
terms  with  various  coefficients,  to  attract  constant  n  lines  to  the  body. 
The  coefficients  have  been  calculated  for  a  flat  plate  to  ensure  sufficient 
constant  n  lines  in  the  boundary  layer.  Those  P  and  Q  functions  were  pro¬ 
grammed  for  airfoil  calculations  in  viscous  flow,  and  are  less  important 
for  potential  flow  calculation  considered  here.  Consequently,  slight 
attraction  of  constant  n  lines  to  the  airfoil  is  used  for  free  stream 
airfoil  calculations  and  Q  is  set  to  zero  in  all  other  cases. 
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Difference  Equations 

The  governing  differential  equations  are  approximated  by  second 
order  central  difference  expressions.  The  number  of  constant  £  lines 
is  designated  IMAX  and  the  number  of  constant  c  lines  by  JMAX,  the  com¬ 
putational  field  is  (JMAX-2) • (IMAX-1 ).  In  the  computer  program  (Appendix 
A)  that  is  used  the  standard  field  size  is  71  x  44.  The  boundary  mapping 
described  previously  is  such  as  J=1  and  J=JMAX  correspond  to  the  inner 
and  outer  boundaries.  The  inner  contour  (the  body)  is  given  by  71  pairs 
of  (x,y)  data  points,  the  outer  boundary  is  either  given  in  the  same 
manner  or  is  calculated  in  the  program.  An  initial  guess  is  generated 
in  the  program.  Because  the  difference  equations  are  nonlinear,  the  ini¬ 
tial  guess  must  be  within  a  certain  neighborhood  of  the  solution  if  the 
iterative  solution  is  to  converge.  Boundary  input  data  and  initial 
guess  type  for  the  various  regions  are  discussed  in  Chapters  V  and  VI. 
Difference  equations  are  solved  by  a  SOR  (Successive  Over  Relaxation) 
iterative  scheme.  The  acceleration  parameter  of  the  SOR  scheme  was  usual¬ 
ly  set  equal  to  one,  since  the  objective  of  this  study  is  not  optimization 
of  transformation  computer  time. 

Transformation  Properties 

The  transformation  procedure  described  above  holds  several  benefits 
and  properties. 

1.  Numerical  solution  of  partial  differential  equations  may  be 
done  in  the  transformed  plane  and  on  a  fixed  rectangular  grid  with  square 
mesh.  Even  though  the  partial  differential  equation  is  complicated  by 
the  transformation,  the  numerical  scheme  does  not  have  to  compensate  for 
grid  changes. 
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2.  No  interpolation  is  required  to  locate  boundary  points  regard¬ 
less  of  the  physical  field.  This  property  is  of  particular  importance 
for  boundaries  with  strong  curvature  such  as  an  airfoil  leading  edge. 

3.  Even  time  dependent  boundaries  transform  to  a  fixed  rectangular 
grid  with  square  mesh,  but  the  transformation  equations  (23  and  (24)  must 
be  re-solved  at  each  time  step,  which  changes  a,3,y,J  as  well  as  the 
coordinate  system  in  the  physical  plane  at  each  time  step. 

4.  The  coordinate  system  generated  by  equations  (23)  and  (24)  is 
not  necessarily  orthogonal  but  it  is  a  general  solution  for  arbitrary 
spacing  and  capability  of  concentration  of  coordinate  lines.  This  gene¬ 
rality  outweighs  the  lack  of  orthogonality. 

5.  The  partial  differential  equation  type  (elliptic,  parabolic, 
hyperbolic)  remains  invariant  under  the  transformation. 

Applications 

Generation  of  specific  boundary  fitted  coordinated  systems  for  the 
various  boundaries  chosen  to  solve  the  flow  field  for  an  airfoil  in  both 
free  stream  and  ground  effect  are  described  in  Chapter  V.  The  coordinate 
system  for  the  hovercraft  model  is  described  in  Chapter  VI. 

As  a  simple  verification  of  the  computer  program,  two  boundary 
fitted  coordinate  systems  were  generated  for  the  region  between  two  con¬ 
centric  circles  of  one  and  ten  units  diameter.  Equally  spaced  concentric 
circles  were  used  as  an  initial  guess  in  each  case.  The  first  coordinate 
system  was  computed  without  attraction  of  coordinate  lines  (P  =  Q  =  0 
in  equation  (23)  and  (24)),  where  the  expected  and  derived  coordinate 
system  consisted  of  concentric  circles  (Figs  8  and  9).  The  same  region 
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is  solved  the  second  time  with  slight  attrection  to  the  inner  circle  lead- 
ing  edge  (point  (-0.5,0)),  the  attraction  can  easily  be  noticed  in  Fig  10 
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FIGURE  8 


Body-Fitted  Coordinate  System  for  a  Cylinder 
in  a  Circular  Outer  Boundary  (Entire  region). 
Each  fifth  s  and  n  line  are  shown. 
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FIGURE  9.  Body  Fitted  Coordinate  System  for  a  Cylinder  in  a 

Circular  Outer  Boundary.  No  Attraction  of  Coordinates 
(P  =  Q  =  0) 
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FIGURE  10.  Body-Fitted  Coordinate  System  for  a  Cylinder  in  a 
Circular  Outer  Boundary.  Attraction  of  Constant  n 
Coordinates. 
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IV.  STREAM  FUNCTION  SOLUTION 


Any  set  of  partial  differential  equations  may  be  solved  on  the 
boundary  fitted  coordinate  system  by  transforming  the  equations  and  asso¬ 
ciated  boundary  conditions  and  solving  numerically  in  the  transformed 
plane.  In  particular,  the  method  is  applied  here  to  solve  the  stream 
function  equation. 

Differential  Equation  and  Aerodynamic  Coefficients 

Two  dimensional  potential  flow  about  an  airfoil  is  described  by 


Laplace  equation 

for  the  stream  function. 

ill  +  ij; 

xx  ryy 

=  0 

(25) 

The  simplest  boundary  conditions  are  Dirichlet  type  where 

are  specified  on  the  boundaries 

values 

<i<(x,y)  = 

iln(x,y)  on  the  body 

(26) 

'l'(x.y)  = 

*2(x,y)  on  the  outer  boundary 

(27) 

Equation  (25)  and  boundary  conditions  (26)  and  (27)  after  trans¬ 
formed  are: 


'  2B^n 

+  y\p  +oi p  +  x\pr  =  0 

nn  rn  C 

(28) 

= 

*i^*nmin} 

(29) 

^’nmax^ 

*2U*nn«x) 

(30) 
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For  this  case  of  irrotational  flow  the  Bernoulli  equation  holds  for 
the  entire  field,  and,  for  stream  function  that  is  nondimensional ized 
relative  to  the  airfoil  chord  and  free  stream  velocity,  we  find  on  the 
body  surface 

Cp  =  1  '  (31) 


In  our  case  the  free  stream  is  in  positive  x  direction  and  there¬ 
fore  the  integration  of  Cp  to  obtain  lift  and  quarter  chord  moment  co¬ 
efficients  leads  to: 

^max 

--  J  cpV?  (32) 

Snin 


Mr  *  §)  dt 


(33) 


a»B,Y,a,x,J  in  the  above  equations  are  functions  of  partial  deriva¬ 
tives  of  the  physical  coordinates  (x,y)  with  respect  to  transformed  coordi¬ 
nates  (s,n)  (Ref  4,  Appendix  A).  It  should  also  be  noted  that  using  the 
Bernoulli  equation  to  calculate  pressure  distribution  is  a  correct  approach 
for  an  irrotational  flow  field.  Solution  of  a  flow  field  for  a  hovercraft 
is  different  in  nature  and  the  differences  are  discussed  in  Chapter  VI. 
Difference  Equation 

Equation  (28)  is  approximated  by  second-order  difference  expressions 
on  the  transformed  grid.  For  the  chosen  field  size,  71  x  44,  Dirichlet 
boundary  conditions  are  specified  by 
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=  ipi ( 1 , 1 )  on  the  body  (34) 

( I , J )  =  4>2(Ij44)  on  the  outer  boundary  (35) 

An  explicit  scheme  is  used  to  solve  the  stream  function  field.  The 

difference  equation  corresponding  to  equation  (28)  is 


T  =  2(“i.j + 


(37) 


Where  superscript  n  refers  to  the  most  recent  stream  function  values. 
Starting  from  an  initial  guess  for  stream  function  values  over  the  compu¬ 
tational  field,  iteration  with  equation  (36)  converges  to  the  solution  . 

The  pressure  coefficient  equation  (31)  is  evaluated  on  the  body  and, 
as  the  body  maps  to  0=1  line,  i|»  is  approximated  by  second  order  one-sided 
difference  expression.  Integrating  the  pressure  coefficient  into  lift  and 
quarter  chord  moment  "oefficients  (Equations  (32)  and  (33))  is  accomplished 
by  using  Simpson's  one-third  numerical  integration  algorithm.  For  com¬ 
puter  programs  see  Appendix  A. 

Appl ications 

Stream  function  solutions  for  an  airfoil  in  free  stream  and  ground 
effect  are  given  in  Chapter  V  and  stream  function  contour  plot  are  shown 
too.  Modifications  of  the  above  equations  to  include  multi-valued  points 
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(branch  cut)  and  vorticity  to  accommodate  for  the  hovercraft  peripheral 
jet  flow  field  solution  are  introduced  in  Chapter  VI. 

To  check  the  equations  and  computer  programs  that  are  used,  the 
pressure  coefficient  variation  for  a  cylinder  in  uniform  flow  is  computed 
using  the  coordinate  system  generated  as  described  in  the  previous  chapter. 
As  shown  in  Fig  11,  good  correlation  is  achieved  between  the  analytical 
and  the  finite  difference  solutions.  Stream  function  contour  lines  for 
this  case  are  shown  in  Fig  12. 


29 


PRESSURE  CGEF 

-3.20  -2.40  -1.60  -0.80  0.00  0.80  1.60 


0.  00 
X/D 


FIGURE  11.  Pressure  Coefficient  on  a  Cylinder  vs.  Nondimensional ized 
Diameter  in  Free  Stream  Direction 
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V.  AIRFOIL  STREAM  FUNCTION  SOLUTION 


This  chapter  describes  the  numerical  solution  for  airfoil  in  free 
stream  and  ground  effect.  Included  are  details  of  the  generation  of 
boundary  fitted  coordinate  system,  the  determination  of  a  stream  function 
solution  with  the  Kutta  condition  satisfied  at  the  leading  edge,  and  the 
derivation  of  lift  and  quarter-chord  moment  coefficients. 

Airfoil  Type 

The  procedure  used  is  applicable  to  any  arbitrary  airfoil  for  which 
the  flow  remains  attached.  For  simplicity  and  to  permit  a  comparison  with 
results  derived  by  finite  singularity  method  (Chapter  II),  a  thin  airfoil 
was  chosen.  The  finite  singularity  method  was  solved  for  a  flat  plate, 
therefore  for  stream  function  solution  a  symmetric  1%  thickness  ratio  air¬ 
foil  is  used.  The  determination  of  the  grid  points  distribution  around 
the  airfoil  without  interpolation  between  data  points  can  be  achieved  by 
using  a  NACA  four-digit  airfoil  that  can  be  described  by  an  equation  (Ref 
5).  Consequently,  the  chosen  airfoil  can  be  designated  as  NACA  0001. 

Kutta  Condition 

Ordinarily,  the  Kutta  condition  at  the  trailing  edge  of  an  airfoil 
can  be  stated  in  terms  of  either  velocities  or  pressure  differences.  The 
most  convenient  manner  of  stating  Kutta  condition  in  our  case  is  in  terms 
of  the  stream  function.  Specifically,  the  body  stream  function  must  detach 
from  the  airfoil  at  the  trailing  edge  with  the  local  camber  line  slope. 

The  procedure  for  stream  function  solution  with  Dirichlet  boundary  con¬ 
ditions  requires  that  a  value  for  the  body  stream  function  be  specified 
and  an  arbitrary  value  will  not  necessarily  satisfy  the  Kutta  condition. 


Thus,  the  computation  of  a  stream  function  solution  which  does  satis¬ 
fy  the  Kutta  condition  requires  an  iteration  and  interpolation  upon  the 
value  chosen  for  the  body  streamline.  Three  points  are  chosen  from  the 
field  solution;  one  is  the  trailing  edge,  and  the  other  two  are  slightly 
behind  the  trailing  edge  and  located  on  opposite  sides  of  the  extension  of 
trailing  edge  camber  line.  The  physical  coordinates  and  the  stream  function 
values  are  known  in  these  points  that  are  marked  as  1 ,  2  and  3  (Fig  13). 


Trailing  Edge  Camber 
Line  Slope 

Figure  13.  Kutta  Condition  Check  Scheme 

Equation  of  the  plane  passing  through  the  three  points  is 

X  Y  1 

Xj  Y  i  'P  i  1 

=  0 

Xj  Y2  t|>2  1 

X3  Y3  *3  1 


(39) 


Substituting  <|<  =  in  equation  (39)  the  plane  equation  reduces  to 
a  straight  line  equation  passing  through  the  trailing  edge  with  the 
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stream  function  value  constant  along  the  straight  line.  The  slope  of 
the  straight  line  is  given  by 


The  Kutta  condition  is  said  to  be  satisfied  whenever  the  slope  derived  by 
equation  (40)  is  in  the  proximity  of  the  trailing  edge  camber  line  slope 
within  specified  allowable  error. 

Airfoil  in  Free  Stream 

The  airfoil  is  the  inner  boundary  for  generation  of  body  fitted 
coordinate  system.  The  point  collection  describing  the  airfoil  contour 
is  calculated  by  the  equation  for  NACA  0001  airfoil  (Table  3).  The  ori¬ 
gin  for  the  region  is  at  the  mid-chord  Which  is  nondimensionalized  to  unit 
length.  In  order  to  use  the  same  airfoil  point  distribution  for  free 
stream  and  ground  effect  cases,  the  airfoil  is  rotated  clockwise  by  the 
angle  of  attack.  The  boundary  fitted  coordinate  system  has  to  be  re¬ 
generated  for  each  angle  of  attack.  The  "infinity11  outer  boundary  is  a 
circle  of  ten  chord  length  radius,  which  is  described  by  equally  spaced 
points  on  the  circumference.  The  airfoil  points  are  numbered  from  the 
trailing  edge  clockwise  and  the  points  on  the  outer  boundary  are  numbered 
for  point  (10,0)  also  clockwise.  The  outer  boundary  is  generated  by  a 
program  subroutine.  The  initial  guess  is  a  series  of  elliptical  n  lines 
with  an  increasing  interval  in  axial  direction.  A  typical  body  fitted 
coordinate  system  is  shown  in  Fig  14  and  an  enlargement  and  more  detailed 
picture  of  the  airfoil  neighborhood  is  given  in  Fig  15. 

The  stream  function  is  computed  in  the  transformed  plane  from  the 
transformed  Laplace  equation.  Stream  function  values  for  a  uniform  flow 
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.  Body  Fitted  Coordinate  System  for  Airfoil  in  Free 
Stream.  Each  fifth  £  and  n  lines  are  shown. 


FIGURE  15.  Body  Fitted  Coordinate  System  for  Airfoil 
in  Free  Stream,  Partial  Region 
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FIGURE  16.  Stream  Function  .Contour  Plot  for  Airfoil 
10°  Angle  of  Attack  in  Free  Stream. 
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at  zero  angle  of  attack  and  a  guess  for  the  body  stream  function  value 
are  used  as  boundary  conditions.  The  initial  guess  is  also  the  stream 
function  field  for  a  uniform  flow  at  zero  angle  of  attack.  A  second  guess 
for  the  body  stream  function  value  is  used  and  the  stream  function  fields 
are  solved  for  those  two  guesses.  The  extent  to  which  the  solution  cor¬ 
responding  to  each  body  stream  function  value  satisfies  the  Kutta  condition 
is  noted  for  each  case.  Then,  assuming  that  the  discrepancy  in  Kutta  con¬ 
dition  is  a  linear  function  of  body  stream  function  value,  a  new  value  for 
body  stream  function  value  is  determined.  Thus,  an  iterative  procedure 
with  linearly  interpolated  values  for  body  stream  function  is  started. 
Iterations  are  said  to  converge  to  the  final  result  when  the  stream  func¬ 
tion  slope  at  the  trailing  edge  equals  the  airfoil  camber  line  slope  at 
this  point  within  allowable  error.  Contour  stream  function  plots  (Fig  16) 
illustrates  satisfaction  of  Kutta  condition. 

Different  grid  point  distributions  have  been  tried  for  airfoil  rep¬ 
resentation.  The  very  first  distribution  with  equally  spaced  chordwise 
intervals  did  not  follow  the  curvature  of  the  leading  edge  and  produced  a 
sharp  leading  edge.  This  case  was  not  used.  The  leading  edge  was  next 
smoothed  by  distributing  11  out  of  a  total  of  71  points  in  the  0.05  chord 
length  from  the  leading  edge.  The  remaining  parts  were  equally  spaced  be¬ 
yond  this  point.  Thus  the  equation  for  X  location  of  grid  points  for  an 
airfoil  in  zero  angle  of  attack  is: 

X*  =  0.5  -  -  (k-1 )  k  =  1 . 30 

=  0.45  -  (k-30)  k  =  31,... 36  (41) 
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The  same  X  locations  are  used  to  calculate  grid  points  on  the  upper  and 
lower  surfaces.  The  lift  coefficient  derived  for  this  point  distribution 
is  21%  less  than  the  analytical  thin  wing  solution,  because  there  are  not 
enough  points  in  the  proximity  of  the  leading  edge  to  follow  the  strong 
variations.  In  order  to  attract  more  grid  points  toward  the  leading 
edge  (X  =  -0.5)  the  basic  distribution  from  equation  (41)  is  modified  by 
the  following  equation. 

X  =  (X*  -  0.5)n  -  0.5  (42) 

Aerodynamic  coefficients  derived  for  grid  distributions  that  are  obtained 
by  equation  (42)  for  n  =  2  and  n  =  3,  are  given  in  Table  4.  The  results 
in  Table  III  indicate  that  the  lift  coefficient  is  over-estimated  by  10% 

(n  =  3)  and  5%  (n  =  2).  The  over-estimation  is  due  to  relatively  large 
intervals  in  the  mid  chord  vicinity  and  toward  the  trailing  edge. 


TABLE  4.  Calculated  Aerodynamic  Coefficients  for  10° 
Angle  of  Attack  for  Different  Airfoil  Grid  Point  Distribution 


C* 

C 

m 

^body 

Kutta  Error 

Basic 

0.869 

-0.0675 

-0.22983 

0.52  X  10"5 

n  =  2 

1.160 

-0.00188 

-0.22913 

0.66  X  10-7 

n  =  3 

1.217 

-0.00450 

-0.23484 

0.11  x  io-7 

Final 

1.150 

-0.00224 

-0.22816 

0.48  X  10"7 

Thin  Ming 
Theory 

1.096 

0.0 

- 

For  a 

flat  plate 

the  analytical 

solution  shows 

that  ^  40%  of  the 

lift  is  generated  by  the  0.1  chord  length  from  the  leading  edge  (Appendix  B). 


This  suggests  that  40%  of  the  points  should  be  in  this  interval,  and  this 
can  be  achieved  by  using  equation  (42)  for  n  =  2.  Thus,  the  final  point 
distribution  employs  equation  (42)  for  the  29  points  that  fall  within  the 
0.1  chord  length,  and  the  remaining  points  are  distributed  by  10  linearly 
increasing  intervals  between  0.1  and  0.5  chord  length  and  equal  intervals 
from  mid-chord  toward  the  trailing  edge  (Table  4).  The  final  airfoil 
representation  over-estimates  the  thin  wing  theory  by  only  4.5%.  It 
seems  reasonable  to  assume  that  a  finer  grid  will  provide  a  better  corre¬ 
lation. 

The  pressure  coefficient  is  derived  from  the  stream  function  solu¬ 
tion  using  equation  (31).  As  the  trailing  edge  is  a  re-entrant  point,  the 
pressure  coefficient  is  calculated  twice  for  this  point;  once  as  a  point 
on  the  lower  surface,  and  second  as  a  point  on  the  upper  surface,  using 
one  sided  difference  expressions.  Because  of  residual  errors  in  the  cal¬ 
culations,  slightly  different  values  are  obtained  in  violation  of  Kutta 
condition.  The  mismatch  is  resolved  by  interpolating  the  value  on  the 
lower  surface  using  sixth  order  Lagrange  interpolation  and  substituting 
the  same  value  to  the  upper  surface  point.  Fig  17  shows  a  typical  chord- 
wise  variation  of  pressure  coefficient. 

Lift  and  quarter  chord  moment  coefficients  are  obtained  by  inte¬ 
grating  pressure  coefficient  (Equations  (32)  and  (33)).  For  an  ideal 
symmetric  airfoil,  the  quarter  chord  moment  coefficient  is  zero  and  a  non¬ 
zero  calculated  moment  coefficient  is  a  measure  to  the  error  introduced 
in  pressure  coefficient.  The  calculated  coefficients  are  shown  in  Table 
5  and  a  plot  of  lift  coefficient  variation  with  angle  of  attack  is  in 
Fig  18.  The  figure  shows  a  very  reasonable  5%  correlation  between  thin 
wing  theory  and  finite  difference  results. 
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TABLE  5.  Aerodynamic  Coefficients  for  Free  Stream  as  Function 
of  Angle  of  Attack  -  Finite  Difference  Method 


o 

C* 

Cmc/4 

4>  body- 

Kutta  Error 

2 

0.230 

-0.31  X  10'4 

-0.04512 

0.161  X  10‘3 

4 

0.460 

-0.37  X  10'4 

-0.09018 

0.421  X  10'3 

6 

0.696 

-0.17  X  10'2 

-0.13719 

0.220  X  10”3 

8 

0.924 

-0.21  X  10'2 

-0.18277 

0.140  X  10‘4 

10 

1.150 

-0.22  X  10'2 

-0.22817 

0.483  X  10"5 

TABLE  6.  Aerodynamic  Coefficients  in  Ground  Effect  as  Function 
of  Ground  Distance  and  Angle  of  Attack  -  Finite  Difference  Method 


a 

G.D. 

Cmc/4 

body 

Kutta  Error 

10 

0.25 

1.379 

-0.355  X  10"1 

-0.110 

0.207  X  10“6 

10 

0.50 

1 .123 

-0.825  X  10'2 

-0.133 

0.222  X  10'7 

10 

0.75 

1.190 

-0.415  X  10'2 

-0.152 

0.195  X  10-6 

10 

1.00 

1.162 

-0.107  X  10'2 

-0.165 

0.660  X  10'7 

10 

1.50 

1.142 

-0.483  X  10'3 

-0.185 

0.366  X  10~6 

4 

0.25 

0.644 

-0.199  X  10"1 

-0.0475 

0.237  X  10-6 

4 

0.50 

0.539 

-0.825  X  10'2 

-0.0565 

0.242  X  10'7 

4 

0.75 

0.500 

-0.366  X  10'2 

-0.0622 

0.386  X  10'9 

4 

1 .00 

0.479 

-0.212  X  10"2 

-0.0672 

0.274  X  TO'7 

4 

1.50 

0.466 

-0.147  X  10'2 

-0.0746 

0.203  X  10'6 
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Airfoil  in  Ground  Effect 

For  the  ground  effect  solution  the  same  airfoil  type  and  airfoil 
body  point  distribution  are  used.  The  outer  boundary  is  of  rectangular 
shape  in  which  the  base  is  the  ground  and  the  top  and  sides  are  at 
"infinity"  (ten  chord  units).  Half  of  the  outer  boundary  points  (36)  are 
on  the  ground  base  with  linearly  increasing  intervals  from  the  mid-base 
toward  the  sides  (Table  7).  The  point  at  the  outer  boundary  at  the  same 
6  as  the  trailing  edge  is  on  the  base  about  20%  of  the  distance  from  the 
mid-base  toward  the  right  side.  This  location  is  picked  to  get  a  physical 
coordinate  system  that  remains  nonsingular  around  the  trailing  edge  (Figs 
19  and  20).  The  initial  guess  for  the  transformation  is  equally  spaced 
intervals  in  each  direction  with  special  attention  to  the  region  between 
the  airfoil  and  the  ground  where  each  guess  has  a  gradually  decreasing 
slope  changing  from  angle  of  attack  to  zero. 

The  stream  function  equation  is  solved  by  the  same  Kutta  condition 

iteration  procedure,  as  before.  The  outer  boundary  and  initial  guess 

values  are  of  a  uniform  flow  at  zero  angle  of  attack.  The  C.  and  C 

t  mc/4 

variations  with  ground  distance  for  4  and  10  degrees  angle  of  attack  are 
given  in  Table  6  and  in  Figures  21  and  22  where  they  are  compared  against 
finite  singularity  results.  It  is  noted  that  a  much  better  correlation 
is  achieved  for  4  degrees  angle  of  attack  than  for  10  degrees.  The  mis¬ 
match  can  be  explained  by  noting  that  finite  singularity  method  is  based 
on  thin  wing  theory  that  is  limited  to  small  angles  of  attack  and  hence 
to  small  lift  coefficients.  Another  point  of  comparison  between  finite 
singularity  method  and  finite  difference  stream  function  solution  arises 
from  comparing  the  variation  of  ground  velocity  (Fig  23)  where  good  cor¬ 
relation  is  obtained. 


45 


Summa  ry 

In  this  chapter  airfoil  flow  field  solution  is  obtained  for  the 
stream  function  equation  for  potential  flow.  Here,  the  stream  function 
equation  is  solved  by  finite  difference  techniques  on  a  body  fitted 
coordinate  system.  A  linear  interpolation  scheme  is  developed  to  satisfy 
the  Kutta  flow  condition  at  the  trailing  edge.  A  criterion  for  grid 
point  distribution  over  the  airfoil  is  established.  To  resolve  the 
strong  leading  edge  gradients,  40%  of  the  grid  point  must  be  distributed 
over  0.1  chord  length  about  the  leading  edge.  Variation  of  and  Cm 

a .  c 

with  angle  of  attack  and  ground  distance  agree  well  to  the  variation  de¬ 
rived  from  finite  singularity  method  (Chapter  II). 


46 


m 


Body-Fitted  Coordinate  System  for  Airfoil  in  Ground 
Effect,  Partial  Region 
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GROUND-DISTRNCE/CHQRD 

Lift  Coefficient  vs.  Ground  Distance. Angle  of  Attack 
as  a  Farameter  (Finite  Difference  Method). 

—  -  Finite  Singularity  Method 

£  -  Finite  Difference  Method  -  4° 

^  -  Finite  Difference  Method  -  10° 


GROUND-O I  ST RNCE/CHORD 

Quarter-Chord  Moment  Coefficient  vs.  Ground  Distance. 
Angle  of  Attack  as  a  Parameter  (Finite  Difference  Method) 

—  -  Finite  Singularity  Method 

0  -  Finite  Difference  Method  -  10° 

£  -  Finite  Difference  Method  -  4° 


o 

OJ 


x/c 


FIGURE  23.  Ground  Velocity  vs.  Nondimensional ized  Length 
in  Free  Stream  Direction. 

(4°-Angle  of  Attack,  0.5-Ground  Distance) 

-  Finite  Singularity  Method 
A  Finite  Difference  Method 
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VI.  HOVERCRAFJ  STREAM  FUNCTION  SOLUTION 

The  hovercraft  discussion  is  limited  to  the  fundamental  case  of 
symmetric  hovering,  although  asymetric  configurations  are  also  feasible. 
The  putpose  is  to  emphasize  the  difference  between  this  case  and  airfoil 
stream  function  solution  and  suggest  appropriate  solutions  to  get  a  well 
defined  jet. 

Hovercraft  Model 

The  peripheral  jet  type  hovercraft  for  this  2-D  model  is  a  rectan¬ 
gle  .  of  0.9  X  0.1  units  with  two  0.1  diameters  semi-circular  ends  (Fig 
24).  The  two  peripheral  jets  are  0.1  units  width  and  are  located  sym¬ 
metrically  near  the  rectangle  ends.  The  body  point  distribution  for 
hovercraft  numerical  definition  assures  enough  points  to  resolve  both 
the  two  semi-circular  ends  and  the  peripheral  jets  in  particular  (Table 
8).  The  outer  boundary  .is  a  rectangle  in  which  the  base  is  the  ground. 


branch-cut 


Figure  24.  Hovercraft  Model 
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Figure  25.  Body  Fitted  Coordinate  System  for 
Hovercraft  in  Ground  Effect 


TABLE  8.  Numerical  Data  for  a  Hovercraft  Model 


TABLE  9.  Numerical  Data  for  Outer  Boundary  for  Hovercraft 
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The  point  collection  for  outer  boundary  definition  assures  good  reso¬ 
lution  for  the  jet  outlets  (Table  9).  To  obtain  a  reasonable-looking 
boundary  fitted  coordinate  system  in  the  peripheral  jet  region  with  the 
limited  size  arrays,  the  outer  boundary  sides  and  top  are  as  shown  (Fig 
25).  The  initial  guess  is  rectangles  with  equally  spaced  intervals  in 
each  direction. 

Multi-valued  Stream  Function 

Peripheral  jets  introduce  the  same  problem  for  a  stream  function 
solution  as  a  fluid  source,  for  which  integration  of  transverse  velocities 
along  any  closed  path  around  the  source  results  in  two  stream  function 
values  at  the  same  point.  The  difference  between  the  two  stream  function 
values  is  proportional  to  the  source  strength,  which  will  be  referred  as 
Aij>.  The  stream  function  solution  can  be  made  single-valued  by  choosing 
a  branch  cut  along  an  arbitrary  line  connecting  the  inner  and  outer  boun¬ 
daries.  Stream  lines  values  change  by  an  amount  Aip  as  they  cross  the 
branch  cut.  Explicit  difference  equation  (36)  remains  valid  for  irrota- 
tional  flow  for  every  point  in  the  field  except  those  which  form  differ¬ 
ences  across  the  branch  cut.  For  the  points  on  the  branch  cut  the  differ¬ 
ence  equation  is 


*i,j =  .j-r'n+i  ,j-i  > 


The  explicit  difference  expression  (Equations  (36)  and  (43))  are 
solved  for  the  computational  field  with  Dirichlet  type  boundary  conditions. 
The  jets  are  assumed  to  be  uniform  flow  at  the  nozzles  and  at  the  outlets. 
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Figure  26.  Stream  Function  Contour  Plot  for 

Hovercraft,  Potential  Flow  Solution 
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with  an  assumed  velocity  at  the  nozzles.  Stream  function  values  are  cal¬ 
culated  by  integrating  the  velocities  where  the  ground  and  hovercraft 
base  between  the  jets  are  assigned  to  be  zero.  A  contour  plot  for  the 
stream  function  is  shown  in  Fig  26,  from  which  it  may  be  concluded  that 
solving  an  irrotational  flow  does  not  provide  a  solution  with  a  well- 
defined  jet  even  though  the  stream  function  is  single  valued. 

Vorticity 

In  order  to  obtain  a  better  defined  jet,  vorticity  is  introduced 
in  two  shear  layers  on  the  jet  boundaries.  Vorticity  in  the  flow  field 
changes  the  governing  equation  to  Poisson  type. 

=  -u)  (44) 

For  simplicity  the  vorticity  will  be  calculated  as  function  of  the  stream 
function  by  assuming  an  arbitrary  velocity  profile  across  the  jet  as 
shown  in  Fig  27.  The  chosen  velocity  profile  is  uniform  velocity  across 
the  center  semi-width  tapering  linearly  toward  the  boundaries. 


Figure  27.  Jet  Velocity  Profile 
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The  total  change  of  stream  function  value  across  the  jet  is 

♦*  =  3  UL/4  (45) 

Stream  function  variation  as  function  of  Y. 


ii  =  2  UL(Y/L)2 

UL/8  +  UL ( Y/L  -  1/4) 

3  UL/4  -  2  UL ( Y/L -1  )2 


0  < Y/L  <1/4 
1/4  < Y/L  <3/4 
3/4  <Y/L  <  1 


(46) 


Vorticity  is  calculated  from  equation  (44). 


u  =  -V2tJj=  -4U/L 
0 

+4U/L 


0  <Y/L  <1/4 
1/4  <Y/L  <3/4 
3/4  <Y/L  <1 


(47) 


Equation  (47)  can  be  written  in  a  more  convenient  form  as: 


u>  =  -v2^  =  -4U/L 
0 

+4U/L 


0  <ij>/i{j*  <1/6 

1/6  <4>U*  <5/6 
5/6  <ip/ip*  <1 


(48) 


Equation  (48)  expresses  vorticity  as  a  function  of  stream  function 
only  and  thus  is  a  very  convenient  form.  The  explicit  difference  ex¬ 
pression  corresponding  to  equation  (44)  is: 


i|/.  a  ~  (R.H.S.  of  equation  43)  +  [J.  .2/2(ct 
»  t  J  I  9  J 


(49) 


Equation  (49),  which  includes  vorticity  and  a  branch  cut,  together 
with  Dirichlet  boundary  condition  as  described  above,  produced  a  well  de¬ 
fined  jet  (Fig  28). 
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Figure  28.  Stream  Function  Contour  Plot  for  Hovercraft, 
Vorticity  in  Jet  Edges 


i 
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Pressure  Derivation 

Cushion  pressure  is  a  very  important  parameter  for  hovercrafts. 
Exponential  theory,  the  commonly  used  analytical  method  for  cushion 
pressure  calculations  (Ref  3),  assumes  that  the  jet  between  the  nozzle 
and  the  ground  forms  a  section  of  a  torrus.  All  the  other  analytical 
methods  have  to  assume  the  jet  geometry  as  well.  Actual  jet  geometry 
variation  from  the  analytical  model  result  in  inaccuracy  in  cushion 
pressure  calculated. 

A  finite  difference  method,  rather  than  assuming  jet  geometry, 
solves  for  it  by  using  the  stream  function  solution.  The  pressure  field 
can  be  obtained  by  the  Poisson  equation  (Ref  6). 

^  (5°) 

Equation  (50)  is  usually  solved  with  Neuman  type  boundary  condition  that 
can  be  determined  by  Euler  equations  on  the  boundary.  This  suggested 
approach  for  pressure  field  calculation  has  not  been  tried  as  part  of  this 
study,  and  is  left  for  further  investigation. 

Summary 

In  this  chapter  the  2-D  flow  field  for  an  air  cushion  vehicle  with 
peripheral  jet  is  obtained  from  the  stream  function-vorticity,  Poisson 
equations.  As  this  case  introduces  a  flow  source  within  the  computational 
field  and,  thus,  raises  a  question  of  a  multivalued  stream  function,  a 
branch  cut  is  introduced.  Solution  of  an  irrotational  flow  field  shows 
that  the  jet  spreads  widely.  In  order  to  obtain  a  better  defined  jet, 
vorticity  is  derived  from  an  arbitrary  assumed  velocity  distribution  at 
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the  jet  outlets  and  assumed  to  be  constant  alonq  stream  lines.  The  in¬ 
vestigation  must  be  continued  to  determine  the  influence  of  the  above 
assumptions  on  the  pressure  field. 
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VII.  CONCLUSIONS 


Stream  function-vortici ty  equations  expressed  in  numerical  finite 
difference  formulation  can  be  used  on  a  body  fitted  coordinate  system 
to  obtain  flow  field  for  both  an  airfoil  varying  degrees  of  ground  effect 
and  air  cushioned  vehicle  with  a  peripheral  jet.  The  airfoil  flow  field 
is  solved  for  the  potential  flow  case.  For  lift  and  quarter-chord  moment 
coefficients  at  various  angles  of  attack  and  ground  distances,  good  agree¬ 
ment  is  observed  between  independent  solutions  with  a  finite  difference 
method  and  a  finite  singularity  method.  A  criterion  for  grid  points 
distribution  about  the  leading  edge  is  set  up  to  resolve  the  strong  grad¬ 
ients  in  this  region  properly.  The  stream  function- vorticity  equation  is 
solved  with  branch  cut  needed  for  fluid  sources  in  the  hovercraft.  To 
obtain  a  well  defined  jet,  vorticity  in  the  edges  of  jets  is  essential 
and  can  be  simply  modeled  with  some  success. 

Several  aspects  remain  to  be  investigated  to  complete  the  feasibi¬ 
lity  determination.  The  first  major  aspect  is  to  develop  a  stable  for¬ 
mulation  to  obtain  the  pressure  field  from  the  stream  function  field  by 
solving  the  Poisson  equation  with  appropriate  boundary  conditions.  This 
appears  straightforward.  Also,  pressure  field  sensitivity  to  both  assumed 
jet  vorticity  distribution  and  vorticity  transport,  must  be  evaluated. 

The  results  of  this  work  indicate  that  a  2-D  stream  function- 
vorticity  solution  to  an  air  cushioned  airfoil  appears  feasible.  If  re¬ 
sults  remain  positive  after  last  questions  are  answered,  the  result  will 
be  a  powerful  analytical  tool  for  flow  field  analysis. 
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APPENDIX  A 


COMPUTER  PROGRAMS  DESCRIPTION 


Appendix  A  contains  short  descriptions  on  the  computer  programs 
that  are  used.  Only  the  five  main  programs  are  discussed;  other  short 
programs  that  are  used  are  not  of  special  interest.  Four  out  of  those 
five  programs  were  basically  written  by  the  AF  Flight  Dynamics  Laboratory, 
Wright-Patterson  AFB,  but  were  changed  to  accommodate  for  the  special 
needs  of  this  study.  The  main  objective  of  each  program  is: 


1. 

Program  CORDMN 

-  Generation  of  boundary  fitted  coordinate 

system. 

2. 

Program  ABGMN 

-  Derivation  of  transformation  parameters 

3. 

Program  PSIMN 

-  Stream  function  field  solution. 

4. 

Program  UVP 

-  Derivation  of  velocity  components  and 

pressure  fields. 

5. 

Program  AERDCO 

-  Derivation  of  aerodynamic  coefficients. 

The  structure  of  the  programs  and  some  input  and  output  information  is  des¬ 
cribed  for  each  program  separately. 

Program  CORDMN 

This  program  generates  2-D  body  fitted  coordinate  system  for  any 
arbitrary  number  of  bodies  in  the  region.  Input  data  consists  of  inner 
boundary  and  optionally  the  outer  boundary  (if  not  generated  by  the  pro¬ 
gram),  max  number  of  grid  points  in  each  direction  in  the  transformed 
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plane,  number  of  bodies  in  the  field,  max  number  of  iterations,  acceler¬ 
ation  parameters,  allowable  X  and  Y  errors,  various  header  lines,  print 
and  plot  control  parameters.  Optionally  the  field  generated  by  a  previous 
run  may  be  read  in  instead  of  the  boundaries.  Output  data  consists  of  X 
and  Y  arrays  that  may  be  printed  out  and/or  written  on  disk  (tape  10). 

Subroutine  BNDRY:  Produces  an  outer  boundary  if  not  given  by  input 

data. 

Subroutine  GUESSA:  Produces  an  initial  guess  for  the  entire  compu¬ 

tational  field. 

Subroutine  TRANS:  SOR  iteration  routine  for  field  solution  (X,Y). 

Subroutine  RHS  or  RHSBL:  Assigns  coefficient  for  exponential  series 

that  are  the  right  hand  side  of  transformation 
equations.  Subroutine  RHS  is  an  exponential  de¬ 
cay  right  hand  side  where  subroutine  RHSBL  co¬ 
efficients  are  derived  from  Blasius  flat  plate 
boundary  layer  solution. 

Function  ERROR:  Calculates  largest  field  error  between  two  con¬ 

secutive  iteration  steps. 

Subroutine  SLEQ:  Simultaneous  solution  of  linear  algebraic  equation 

used  to  derive  RHSBL  coefficients. 

Subroutine  MAXMIN:  Calculates  the  maximum  and  minimum  values  of  a 

two  dimensional  array.  Used  for  plots  scaling. 
Subroutine  REED:  Reads  data  from  disks. 

Subroutine  WRITDK:  Writes  on  disk  and  prints  output  data. 

Subroutine  CORPLT:  Plots  first  guess  and  final  coordinate  system. 
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Program  ABGMN 

This  program  calculates  numerically  (second  order  difference  ex¬ 
pressions)  the  transformation  parameters  for  the  coordinate  system  arrays 
(X,Y)  that  were  produced  by  program  CORDMN.  It  should  be  noted  that  the 
calculations  give  4a,  4g,  4y,  4J,  2a,  2t,  and  not  the  transformation  para¬ 
meters  themselves.  Input  data  consists  of  X  and  Y  arrays  that  are  read 
from  disk  (tape  10),  control  parameter  for  sharp  trailing  edge,  control 
parameter  for  sharp  trailing  edge,  control  parameter  for  body  transfor¬ 
mation,  parameters  change,  and  various  print  parameters.  Output  data  is 
written  on  disks  (tape  11  -  4a,  4g,  4y;  tape  12  -  4J;  tape  13  -  2a,  2x). 


Subroutine  ABGd: 

Subroutine  INOUT: 
Subroutine  CHNG: 


Subroutine  TEDGE: 


Subroutine  MAXMIN: 

Subroutine  REED: 
Subroutine  WRITDK: 


Calculates  the  transformation  parameters  by  second 
order  difference  expressions. 

Controls  data  input  and  output. 

Calculates  the  expression  a/J  and  y/J  for  J  =  1 
(on  the  body)  and  stores  the  values  as  3  and  y  on 
the  body  (note  that  B  and  y  on  the  body  are  not 
used  for  stream  function  field  solution). 
Recalculates  transformation  parameters  using  central 
differencing  if  the  body  does  not  have  sharp  trail¬ 
ing  edge. 

Calculates  the  maximum  and  minimum  values  of  a  two 
dimensional  array.  Used  for  plots  scaling. 

Reads  data  from  disks. 

Writes  on  disk  and  prints  output  data. 
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Program  PSIMN 

This  program  produces  a  stream  function  field  solution  by  solving 
explicit  numerical  expressions.  Input  data  consists  of  coordinate  sys¬ 
tem  given  by  X  and  Y  arrays  (tape  10),  transformation  parameters  (Tapes 
11  and  13),  max  number  of  iterations,  convergence  criterion,  angle  of 
attack,  free  stream  velocity,  Kutta  convergence  criterion,  max  number 
of  Kutta  iterations,  stream  function  value  on  the  body,  and  print  con¬ 
trol  parameters.  For  hovercraft  solution  jets  velocity  is  given  as  input. 
Output  data  consists  of  stream  function  field  array  (tape  14)  and  body 
pressure  coefficient  one  dimensional  array  (tape  16). 


Subroutine  BNDRYP: 
Subroutine  GUESSP: 
Subroutine  COEFS: 

Subroutine  LAPSOR: 
Subroutine  PRES: 
Function  ERROR: 

Subroutine  REED: 
Subroutine  WRITDK: 
Subroutine  KUTTA: 

Subroutine  OMEGA: 


Produces  outer  boundary  condition. 

Produces  stream  functional  initial  guess. 

Produces  new  coefficients  needed  by  stream  function 
numerical  expression. 

Produces  field  solution. 

Calculates  pressure  coefficient  by  Bernoulli  equation. 
Calculates  largest  field  error  between  two  consecutive 
iteration  steps. 

Reads  data  from  disks. 

Writes  on  disk  and  prints  output  data. 

Calculates  body  stream  function  slope  at  the  trailing 
edge  to  compare  against  camber  slope  (airfoil  only). 
Calculates  vorticity  values  a  function  of  stream  func¬ 
tion  values. 
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Program  UVP 

This  program  produces  velocity  components  and  invicid  pressure 
fields.  Input  data  consists  of  coordinate  systems  given  by  X  and  Y  arrays 
(tape  10),  stream  function  field  array  (tape  14),  angle  of  attack,  velo¬ 
city  at  infinite  to  reference  velocity  ratio,  control  parameter  to  choose 
between  2-D  and  axi symmetric  case  and  print  control  parameter.  Output 
data  consists  of  u  and  v  velocity  components  and  pressure  fields  (tape  20). 


Subroutine  VEL: 

Subroutine  CPINV: 

Subroutine  REED: 
Subroutine  WRITDK: 


Produces  velocity  component  field  by  numerical  dif¬ 
ferentiation  of  stream  function  field. 

Produces  pressure  field  for  invicid  flow  by  Bernoulli's 
equation. 

Reads  data  from  disks. 

Writes  on  disk  and  prints  output  data. 


Program  AEROCO 

This  program  calculates  lift  and  moment  coefficients  for  an  airfoil 
in  invicid  flow  field  by  Bernoulli  equation.  Input  data  consist  of  coordi¬ 
nate  system  given  by  X  and  Y  arrays  (tape  10),  pressure  coefficient  on  the 
body  (tape  16),  and  angle  of  attack.  Output  data  consist  of  pressure  co¬ 
efficient  printout  and  lift  and  moment  coefficients. 


Subroutine  CLCM: 


Subroutine  LAGINT: 


Subroutine  WRIPR: 


Calculates  lift  and  moment  coefficients  by  inte¬ 
grating  numerically  pressure  coefficient  by  Simpson’s 
one  third  rule. 

Corrects  pressure  coefficient  at  the  trailing  edge 
by  using  Lagrange  interpolation  formula. 

Writes  on  disk  and  prints  single  dimensional  array 
output  data. 
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APPENDIX  B 


Chord-Wise  Cp  Variation  for  Flat  Plate 

Appendix  B  contains  an  analytical  analysis  for  chord-wise  Cp  vari¬ 
ation  for  a  flat  plate.  The  analysis  is  done  using  thin  wing  theory 
vortices  distribution  for  a  flat  plate  (Ref  1).  The  form  of  the  derived 
results  is  percentage  of  lift  coefficient  vs.  non-dimensional ized  chord 
length  from  the  leading-edge.  The  results  are  used  as  a  measure  for 
point  distribution  on  the  airfoil  that  will  pick  up  the  high  gradients 
around  the  leading  edge  without  amplifying  Cp  values  on  the  rear  section 
of  the  airfoil . 

Lift  coefficient  is  calculated  by  integration  of  vorticity  distri¬ 
bution,  which  leads  to  the  following: 

e  -  arccos  (2  x/c-1 )  (B-l ) 

%U  =  100  [1  -  (e-si no  )/n  ]  (B-2) 

the  results  from  B-l  and  B-2  are  given  in  Table  B-l  and  Figure  B-l.  Re¬ 
sults  show  that  40%  of  the  lift  for  a  flat  plate  is  generated  by  I  Oil  of 
chord  from  the  leading  edge,  and  only  18%  of  the  lift  is  generated  by  the 
semi-chord  between  mid-chord  and  trailing  edge.  Conclusion  from  the  lift 
contribution  by  various  parts  of  the  airfoil  is  that  very  rapid  pressure 
changes  and,  hence,  velocitv  changes  occur  in  the  vicinitv  of  the  leadina 
edge.  In  the  finite  difference  numerical  field  solution  where  the  matrices 
sizes  are  tried  to  be  kept  reasonable,  more  points  should  be  put  near  the 
leading  edge  to  pick  up  those  changes.  The  chosen  point  distribution, 
that  show  to  give  good  results,  was  40%  of  the  points  within  10%  from  the 
leading  edge  and  only  27%  of  the  point  on  the  rear  semi-chord. 
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Table  B-l.  Percent  cf  C;  vs.  Nondimensional  Chord 
Measured  from  Leading  Edge 


x/c 

%  C* 

•  x/c 

of  c 

*  4 

0. 

0. 

0.10 

39.6 

0.01 

12.7 

0.20 

55.0 

0.02 

17.9 

0.30 

66.1 

0.03 

21 .9 

0.40 

74.8 

0.04 

25.3 

0.50 

81 .8 

0.05 

28.2 

0.60 

87.6 

0.06 

30.9 

0.70 

92.3 

0.07 

33.3 

0.80 

•  95.9 

0.08 

35.5 

0.90 

98.6 

0.09 

37.6 

1.00 

Figure  B-l.  Graphic  Presentation  of  Table  B-l 
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